subroutine inout
  use rbf_type
  use variable
  use random_number2
  use ranftest


  implicit none
  logical(lgt) antigua,iniciar
  integer(i4b) i
  real(dp) ::  S,p0

  !-- Parameters

   pi=acos(-1.d0)

  !--- input file ---

  open(111,file='drop.in',status='old')
  read(111,nml = input)
  close(111)

  ! Checking for previous 
  inquire( file = 'drop.out', exist = antigua )
  inquire( file = 'rand.bin', exist = iniciar )


  antig:if ( .not. antigua ) then

     open(500,file='results.bin',form='binary',status='new',action='write')
     open(501,file='scalar.dat',status='new')
     open(502,file='biaxial.dat',status='new')
     open(503,file='energy.out',status='new')
     open(504,file='energy_KT.out',status='new')

     if(.not.iniciar)then
        open(11,file='rand.bin',form='binary',status='new')
        call rantest(iseed)
        close(11)
     else
        call restart_rand(iseed)
     endif


  ! Elastic constants and parameters     
  lhat(1) = 1.d0
  LL      = ll1*lhat          ! Dimensions of the other constants
  radius   = beta*0.5d0*nematic_l   ! nm

  ! For bulk relaxation only
  Uinitial = ULdG           ! For equilibration and avoid frozen configurations
  Sinitial = S( Uinitial )
  Sbulk    = S( ULdG )            ! Scalar Order parameter 

  ALdG   = abs( LL(1) )/( nematic_l*1.d-9 )**2          ! Landau de Gennes constant (j/m^3)
  Ws     = abs( LL(1) )/( extrap_ls*1.d-9 )             ! Surface Energy Surfactant (j/m^2)
  Ww1    = abs( LL(1) )/( extrap_lw1*1.d-9 )            ! Surface Energy Water 1     (j/m^2)
  Ww2    = abs( LL(1) )/( extrap_lw2*1.d-9 )            ! Surface Energy Water 2     (j/m^2)

  gamma_s   = nematic_l/extrap_ls                  ! Dimensionless extrapolation lenght Surfactant
  gamma_w1  = nematic_l/extrap_lw1                 ! Dimensionless extrapolation length Water 1
  gamma_w2  = nematic_l/extrap_lw2                 ! Dimensionless extrapolation lenght Water 2
  gamma_w2 = 0.d0                                  ! Following some papers

! Ossen-Frank elastic constants
  kks(1)  = ( 2.d0*LL(1) + LL(2) + LL(4) )*Sbulk**2 - ( 2.d0*LL(3)*Sbulk**3 )/3.d0
  kks(2)  = ( -2.d0*Sbulk**2*( -3.d0*LL(1) + LL(3)*Sbulk ) )/3.d0
  kks(3)  = ( 2.d0*LL(1) + LL(2) + LL(4) )*Sbulk**2 + ( 4.d0*LL(3)*Sbulk**3 )/3.d0
  kks(4)  = LL(4)*Sbulk**2
! Chirality term


  q0 = lhat(5)*ll1*Sbulk**2/(2.d0*kks(2))


  ! Check on Kij inequalities
  if ( kks(1) < 0 ) print *,'-->--> Violation of Ericksen Inequality: K11 > 0 <--<--'
  if ( kks(2) < 0 ) print *,'-->--> Violation of Ericksen Inequality: K22 > 0 <--<--'
  if ( kks(3) < 0 ) print *,'-->--> Violation of Ericksen Inequality: K33 > 0 <--<--'
  if ( kks(2) + kks(4) > 2 )  print *,'-->--> Violation of Ericksen Inequality: 2K11 > k22 + k24 <--<--'
  print *,' '
  write(*,100)'    Splay elastic - K11    (N)  :', kks(1)
  write(*,100)'    Bend elastic  - K22    (N)  :', kks(2)
  write(*,100)'    Twist elastic - K33    (N)  :', kks(3)
  write(*,100)'    Saddle-Splay el. K24   (N)  :', kks(4)
  print *,' '
  write(*,100)'    Radius to Nematic Coherence        :', beta
  write(*,100)'    Radius to extrapolation surfactant :', gamma_s
  write(*,100)'    Radius to extrapolation Water 1    :', gamma_w1
  write(*,100)'    Radius to extrapolation Water 1    :', gamma_w2
  write(*,100)'    Surface Energy Surfactant (j/m^2)  :', Ws
  write(*,100)'    Surface Energy water      (j/m^2)  :', Ww1
  print *,' '
  write(*,110)'    Initial Nematic-Isotropic parameter  :', Uinitial
  write(*,110)'    Initial bulk scalar order parameter  :', Sinitial
  write(*,110)'    Final Nematic-Isotropic parameter    :', ULdG
  write(*,110)'    Final bulk scalar order parameter    :', Sbulk

  write(*,120)'    AldG Landau de Gennes (j/m^3) (Zumer): ',AldG*(1.d0 - ULdG/3.d0)
  write(*,110)'    Nematic Coherence L. sqrt(L1/AldG)   :',nematic_l



  if(Nq0.ne.0.d0)then
     p0 = (4.d0*(radius)/Nq0)
     q0 = 2.d0*pi/(p0/nematic_l)
     lhat(5) = (2.d0*kks(2)*q0/Sbulk**2)
     LL(5)      = lhat(5)
     lhat(5)    = lhat(5)/ll(1)
  
     if(lhat(5).ne.0)then
        write(*,*)
        write(*,110)'    Pitch (nm)  :', p0
        write(*,110)'    q0          :', q0
        write(*,110)'    lhat (5)    :', lhat(5)
        
     endif
! q0 = lhat(5)*ll1*Sbulk**2/(2.d0*kks(2))

  else
     if(lhat(5).ne.0)then
        write(*,110)'    Pitch (nm)  :', nematic_l*2.d0*pi/q0
        write(*,110)'    q0          :', q0
     endif
  endif


print *,' '
100 format (a,es12.4)
110 format (a,f10.4)
120 format (a,es12.4)


    call grid_uniform

    print *,'   Boundary Nodes           :', nn
    print *,'   Random internal points   :', nip
    print *,'   Total Number of Unknowns :', nn + nip 
    write(*,*)'------------------------------------------------------------'
    print *,' '

    ! dimension variables
    nodes = nn + nip 

    call cube
    
    allocate( normal(3*nodes), planar(3*nodes) )
    allocate (tangential(3*nn),tangential_2(3*nn))
    allocate(  a(5,nodes), da(3*nodes,5),q(5,nodes))
    allocate(  chi(nn), ao_s(5,nn),ao_w(5,nn))

    call drop_normal

    call q_anchoring
    call q_initial     ! --> Design tensors for the bulk
    call initial_da
   

    write(500) xm,  a, nee_cube,ao_s,ao_w,chi,inside,surface

    call grid_nano
    if(npatch.ne.0)then
       call nano_normal
       call q_nano 
       
    endif


    vi = (8.d0*ry*rz*beta**3)  /nip
    ai = (8.d0*ry*beta**2)/(nn )

    close(500)
    close(501)
    close(502)
    close(503)
    close(504)

     !--- Output namelist ---
     open(111,file='drop.out',status='unknown')
     write(111,nml = addic)
     close(111)

  else
  

     !--- Output namelist ---
     open(111,file='drop.out',status='unknown')
     read(111,nml = addic)
     close(111)

     open(500,file='results.bin',form='binary',status='old',action='READ')

     
     ! Principal memory
        allocate( xm(3*nodes) )
    
     ! Memory All nodes
        allocate(  a(5,nodes),nee_cube(Nx,Ny,Nz),da(3*nodes,5),q(5,nodes))
        allocate (tangential(3*nn),tangential_2(3*nn))
        allocate( normal(3*nodes), planar(3*nodes) )
        allocate(  chi(nn), ao_s(5,nn),ao_w(5,nn))
        allocate(  inside(Nx,Ny,Nz),surface (Nx,Ny,Nz))

   
        read(500) xm,  a, nee_cube,ao_s,ao_w,chi,inside,surface


        do i = 1, archivo_cop
           read(500) a, energy        
        end do

        call cube
        call drop_normal

        nip = nodes-nn
    
        ddx = xmax/dble(Nx-1)
        ddy = ymax/dble(Ny-1)
        ddz = zmax/dble(Nz-1)

        dddx = 0.5d0/ddx 
        dddy = 0.5d0/ddy
        dddz = 0.5d0/ddz

        ddx2 =  1.d0/ddx**2
        ddy2 =  1.d0/ddy**2
        ddz2 =  1.d0/ddz**2

        

        call initial_da
!**************************************
        call grid_nano
        if(npatch.ne.0)then
           call nano_normal
           call q_nano 
       endif

!-------------------------------------
        char_ener     = ALdG*( nematic_l*1.d-9 )**3

        if(Nq0.ne.0.d0)then
           p0 = (4.d0*(radius)/Nq0)
           q0 = 2.d0*pi/(p0/nematic_l)
           lhat(5) = (2.d0*kks(2)*q0/Sbulk**2)
           LL(5)      = lhat(5)
           lhat(5)    = lhat(5)/ll(1)
           
           if(lhat(5).ne.0)then
              write(*,*)
              write(*,110)'    Pitch (nm)  :', p0
              write(*,110)'    q0          :', q0
              write(*,110)'    lhat (5)    :', lhat(5)
              
           endif
           
        else
           if(lhat(5).ne.0)then
              write(*,110)'    Pitch (nm)  :', nematic_l*2.d0*pi/q0
              write(*,110)'    q0          :', q0
           endif
        endif

        vi = (8.d0*ry*rz*beta**3)  /nip
        ai = (8.d0*ry*beta**2)/(nn )


       close(500)

    
    endif antig



end subroutine inout
!-----------------------------------------------------------
!-----------------------------------------------------------
!-----------------------------------------------------------
subroutine output
  use rbf_type
  use variable, only: ALdG,nematic_l,radius,q, a, energy,iseed,addic,plan,nn,ddz,zmax
  use variable, only: archivo_cop,ini_step ,timesteps,step,nodes,ry,rz,nee_cube,nx,ny,nz,pi

  use random_number2
  use ranftest

  implicit none
  integer(i4b) i,ctemp,iii,k,j,count_spin
  Character(len=10) :: namefile
  Character(len=20), Parameter :: sfile_1= 'drop.out.'
  real (dp) :: aux1,aux2,aux3,scalar,scalar_ave,eta_ave,eta,norm,angle,ave_eta,ave_s,ave_azimu_ang,angle_azi,ave_polar_ang,angle_pol
  real(dp),dimension(5) :: q_tmp
  real(dp),dimension(3) :: director
  logical :: arch
  save ctemp


  aux1 = ALdG*( nematic_l*1.d-9 )**3               ! Characteristic Energy
!  aux2 = 4.0*pi*( radius*1.d-9 )**3/3.d0           ! Volume
!  aux2 = (8.d0*ry*rz*(radius*1.d-9)**3)  
!  aux3 = 249.35d0/1000.d0/1000.d0                  ! Molecular weight 
!  aux1 = aux1/aux2*aux3      



  open(500,file='results.bin',form='binary',status='old',access='append',action='WRITE')
  write(500)  a, energy
 
  archivo_cop = archivo_cop + 1


  if(step.eq.timesteps) ini_step = timesteps

  !--- Output namelist ---
  open(111,file='drop.out',status='unknown')
  write(111,nml = addic)
  close(111)
  close(500)
  call rand_save(iseed)
  
!  Save configuration when change temperature

  if (step.eq.0) ctemp = 0
  if (step.gt.0.and. mod(step,timesteps/5).eq.0) ctemp = ctemp + 1
 !-------------------------------------------------------------------------------------
  open(501,file='scalar.dat',status='old',POSITION='append')
  open(502,file='biaxial.dat',status='old',POSITION='append')
  open(503,file='energy.out',status='old',POSITION='append')
  open(504,file='energy_KT.out',status='old',POSITION='append')
!------------------------------------------------------------------
  open(505,file='angle.dat',            status='replace')
  open(506,file='eta_profile.dat',      status='replace')
  open(507,file='s_profile.dat',        status='replace')
  open(508,file='azimuthal_profile.dat',status='replace')
  open(509,file='polar_profile.dat',    status='replace')
  open(510,file='current_scalar.dat',   status='replace')
  open(511,file='current_biaxial.dat',  status='replace')
! -------
  write(503,100)step, energy
  write(504,101)step, energy*aux1/kbt

  scalar_ave = 0.d0
  eta_ave    = 0.d0

  do i=1,nodes
     call inv_transformation( a(:,i), q_tmp (:))
     call autovalores_biax(q_tmp(:), scalar, eta,director(:))
     scalar_ave = scalar_ave  + scalar*1.5d0
     eta_ave    = eta_ave + (eta*3.d0 + scalar*1.5d0)
     write(510,120)i,scalar*1.5d0
     write(511,120)i,(eta*3.d0 + scalar*1.5d0)
  enddo

  write(501,*)step,scalar_ave/dble(nodes)
  write(502,*)step,eta_ave/dble(nodes)

  do iii=1, nz
     ave_eta       = 0.d0 
     ave_s         = 0.d0
     ave_azimu_ang = 0.d0
     ave_polar_ang = 0.d0
     count_spin    = 0

     do j=1, ny
        do k=1, nx
          
           i=nee_cube(k,j,iii)

           call inv_transformation  (a(:,i), q_tmp(:))
           call   autovalores_biax  ( q_tmp(:), scalar,eta, director(:))

           if(abs(director(2)).le.1.d-8)then
              if(director(1).lt.0.d0)then
                 angle_azi = -pi*0.5d0
              else
                 angle_azi = pi*0.5d0
              endif
           else
              angle_azi = atan(director(1)/director(2))
           endif

           if(abs(director(3)).le.1.d-8)then
              if(director(2).lt.0.d0)then
                 angle_pol = -pi*0.5d0
              else
                 angle_pol = pi*0.5d0
              endif
           else
              angle_pol = atan(director(2)/director(3))
           endif

           ave_azimu_ang = ave_azimu_ang + abs(angle_azi)
           ave_polar_ang = ave_polar_ang + abs(angle_pol)
           ave_eta = ave_eta + (eta*3.d0 + scalar*1.5d0)
           ave_s = ave_s + scalar*1.5d0

           count_spin = count_spin + 1
        enddo
     enddo

     write(506,121)(ddz*(iii-1)/(zmax)-0.5 ),ave_eta/dble(count_spin)
     write(507,121)(ddz*(iii-1)/(zmax)-0.5 ),ave_s/dble(count_spin)
     write(508,121)(ddz*(iii-1)/(zmax)-0.5 ),180.d0*(ave_azimu_ang/dble(count_spin))/pi
     write(509,121)(ddz*(iii-1)/(zmax)-0.5 ),180.d0*(ave_polar_ang/dble(count_spin))/pi
 
  enddo


  close(501) ; close(502) ; close(503) ; close(504)
  close(505) ; close(506) ; close(507) ; close(508)
  close(509) ; close(510) ; close(511)

100 format (I8,8(2x,es16.8))
101 format (I8,8(2x,es16.8))
120 format (I8,es16.8)
121 format (2f12.6)

end subroutine output
!-----------------------------------------------------------
!-----------------------------------------------------------
!-----------------------------------------------------------
